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ABSTRACT 

A detailed description of spectral multigrid methods is 
provided. This includes the interpolation and coarse-grid 
operators for both periodic and Dirichlet problems. The spectral 
methods for periodic problems use Fourier series and those for 
Dirichlet problems are based upon Chebyshev polynomials. An 
improved preconditioning for Dirichlet problems is given. 
Numerical examples and practical advice are included. 
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I. INTRODUCTION 


The motivation for applying a pseudospectral discretization 
to elliptic problems is to obtain an highly accurate app "oxi- 
mation with a small number of collocation points. The major 
advantage that this sort of discretization often offers over 
standard finite difference or finite element techniques is 
greatly reduced storage requirements. At the NASA Ames Symposium 
on Multigrid Methods, we proposed a spectral multigrid approach 
to solving the discrete equations which arise . from applying 
pseudospectral . approximations to variable-coefficient, self- 
adjoint elliptic equations [1]. The focus of that preliminary 
report was on problems with periodic boundary conditions. We 
demonstrated that the number of multigrid iterations necessary to 
achieve convergence was independent -Of the size of the problem. 
The tentative results given in [1] for problems with Dirichlet 
boundary conditions were not so satisfactory because the required 
number of multigrid iterations increased with the number of grid 
points. The purpose of this paper is to fill in some of the 
details omitted from [1] because of page constraints, to describe 
an improved version of spectral multigrid for Dirichlet problems, 
and to offer some practical advice for implementation. 


II. SPECTRAL MULTIGRID ON A SIMPLE MODEL PROBLEM 
The fundamentals of spectral multigrid (SMG) are perhaps 
easiest to grasp for the simple model problem 

- = f (1) 
dx 2 
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on (0 , 2ir ) with periodic boundary conditions. We will examine 
this trivially-solvable problem in detail for the benefit of 
those unfamiliar with either spectral or multigrid methods. The 
standard collocation points are 

Xj = ill j = 0,1, . . . ,N-1 . (2) 


Let fj = f ( x j ) and let Uj be the approximation to u(xj). The 
discrete Fourier, coefficients of Uj are 


u 



Ui6 -2TTijp/N 
j=0 3 




The inverse relationship can be written 

N _ i 

u. = ^ l ue^j . (4) 

J N p 

P = " 2 

Thus, a sensible approximation to the left-hand side of Eq. (1) 
at the collocation points is 

N _ i 

^ l p* u e ipx j . (5) 

N p 

P = ‘ 2 

The pseudospectral approximation to Eq. (1) may be represented by 


where 


LU = F, 


( 6 ) 


u - (Uq,U^ , . . . ,Ujg_^ ) , 


(7) 
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( 8 ) 


p — ( ' ^1 • * * * ^N-l ) • 


and 


L = C“ 1 DC . (9) 

The matrix C represents the discrete Fourier transform; its 
elements are 


_ 1 -2irip j/N 
Pj A 


Clearly, 


(C" 1 ). = e 2irijp/N . 

IP 


( 10 ) 


( 11 ) 


The diagonal matrix D represents the second derivative in trans- 
form space: 




( 12 ) 


In Eqs . (5), (10), (11), and (12) the indices p and q have the 

range indicated in Eq. (3); refer to Eq. (2) for the range of j. 

A Richardson's iterative scheme [2] for solving Eq. (6) is 

V «- V + w(F - LV), (13) 

where V is the current approximation to U and uj is a relaxation 
parameter. The eigenfunctions of L are 


- 4 



« a (p> = e 2 " i 3P/ R , 


(14) 


with the corresponding eigenvalues 

X(p) = P 2 • (15) 

The ranges of j and p are the same as above . The index p has a 
natural interpretation as the frequency of the eigenfunction. 

The error at any stage of the iterative process is V - U; it 
can be resolved into an expansion in the eigenvectors of L. Each 
iteration reduces the p'th error component to v(X^) times its 
previous value, where 

v ( X) = 1 - u,X . (16) 

The optimal choice of & results from minimizing 


I v( X) I for Xe [x . , X 1 , where 

*• m i n ms v J 


X • = 1 and X__„ = N 2 /4 ■ 

min max ' 


min' max- 

(One need not worry about the p = 0 eigenfunction since it corre- 
sponds to the mean level of the solution, which is at one's 
disposal for this problem. ) The optimal relaxation parameter for 
this single-grid procedure is 


J SG " 1 + X . • 

max mm 


(17) 


It produces the spectral radius 


p_^ *max *min . 

SG = x nr~ 

max min 


(18) 


Unfortunately, p__ =* 1 - 8/N 2 , which implies that 0(N 2 ) iter- 

SG 

ations are required to achieve convergence. 

This slow convergence is the outcome of balancing the 
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damping of the lowest-frequency eigenfunction with that of the 
highest-frequency one in the minimax problem described after Eq. 
(16). The multigrid approach takes advantage of the fact that 
the low-frequency modes (|p|‘ < N/4) can be represented just as 
well on coarser grids. It settles for balancing the middle- 
frequency eigenfunction (|p| = N/4) with the highest-frequency 

one ( | p| = N/2), and hence damps effectively only those modes 
which cannot be resolved on coarser grids. In Eqs. (17) and 
(18), ^ m i n replaced with = X(N/4). The optimal relax- 

ation parameter in this context is 


i l. 


MG 


J =r“X"TT 

max mxd 


The multigrid smoothing factor 


(19) 


^MG 


*max “ ^mid 

1 + X " 

max mid 


( 20 ) 


measures the damping rate of the high-frequency modes. In this 
example u.._ = 0.60, independent of N. The price of this 
effective damping of the high-frequency errors is that the low- 
frequency errors are hardly damped at all. However, on a grid 
with N/2 collocation points, the modes for I pi e [n/8, n/ 4] are 
now the high-frequency ones. They get damped on this grid. 
Still coarser grids can be used until relaxations are so cheap 
that one can afford to damp all the remaining modes, or even to 
solve the discrete equations exactly. 

The spectral multigrid approach requires the use of a 
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sequence of grids (or levels). Denote these levels by the index 
k where k = 2,3,...,K. Level k consists of collocation 
points, where = 2 k . Equations (2) to (12) apply on each level 
with replacing N. Throughout this paper the symbol k will be 
used solely to denote a multigrid level. 

On each level we must define a discrete problem, a relax- 
ation scheme, and interpolation operators. The discrete problems 
will be denoted by 


L^V* = 



( 21 ) 


On the finest level, L k = L, F k = F, and V k = U, the solution to 
Eq. (6). The relaxation scheme will be confined here to be 
Richardson iteration 


v k «■ v k + w k (F k - L k v k ) 


( 22 ) 


where v k is the approximation to V k and w k is the relaxation 
parameter. The interpolation operator R k represents the fine-to- 
coarse restriction of residuals from level k to level k - 1: 


F k- 1 = R k (F k - L k v k ). (23) 

The interpolation operator P k represents the prolongation of 
corrections from level k-1 to level k: 


v k + p k v k-l 


(24) 
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Appropriate choices for the coarse-grid operators and the 
interpolation operators are discussed in the next section. Many 
choices are possible for the scheduling algorithm which controls 
the transfer between grids. We will return to this issue in the 
section on numerical examples. 

III. INTERPOLATION AND COARSE GRID OPERATORS 
We will focus on the one-dimensional problem 



on either (0,2it), as in the periodic case or on (-1,1), as in the 
Dirichlet case. We will occasionally refer to equations from 
[1], denoting them by the prefix I, e.g. Eq. (1.5). 

Fourier Series 

The natural interpolation operators represent trigonometric 
interpolation. They were defined in [1] by Eqs. (1.31) and 
(1.32). Useful explicit representations of the restriction and 
prolongation operators (with the superscript k suppressed) are 



1 

N 



e 2uiq(2 j— £ )/N 


(26) 


and 



2 

n 


(27) 




e 2niq(j - 2£)/N 


These summations may be performed in closed form to yield 




(28) 


and 


j* 


2 S 


~ 2 j 


(29) 


where 


^ - 1 r = 0 (mod N) 

sin (^) cot (^)- cos(Z^) otherwise 


(30) 


In analyzing the coarse-grid operator, Eqs . (26) and (27) are 
more useful than Eqs. (28) -(30). Moreover, as noted in Hi], the 
interpolation can be implemented efficiently by Fast Fourier 
Transforms rather than by using Eqs. (28)- (30) in matrix-vector 
multiplications. By the way, the definition of C given here in 
Eq. (10) differs slightly (by a factor /N ) from the definition 
used in [1]. Note that except for a factor of 2, P and R are 
adjoint . 

The pseudospectral evaluation of the left-hand side of Eq. 
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(25) can be expressed as 


MAM U, (31) 

where 

M = C -1 DC , (32) 

A. = a (x . ) 6 . , (33) 

jl 3 j,l 

and in a slight change of notation the diagonal matrix D which 
represents the first derivative in wavenumber space is given by 



(34) 


The reason for setting Dpp = 0 for p = - N/2 is given in [1]. 
Equation (31) costs only 0(N Jtn N) operations to evaluate when 
the Fast Fourier Transform is used. A simple, efficient, and 
effective choice for the coarse-grid operator is 


L^-l = j^k-l^k-l^k-l f 


(35) 


where A k is the diagonal matrix given by 

A* = a k (x.)S. „ , (36) 

jl j 3 , % 

with a K (Xj) = a(Xj), and for k = 3,4,...,K, 
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( 37 ) 


~k-l „k~k 
a = R a 


In other words, the variable coefficient to be used on the 
coarser grid k-1 is a filtered version of the coefficient on 
level k. Otherwise, the coarse-grid operators are the natural 
pseudospectral approximations on those levels. 

As has been stressed especially by Nicolaides [3], Hackbusch 
C4], and Wesseling [5], it seems desirable to use 

L k-1 = R k L k P k , (38) 

with R k the adjoint of P k . The choice made above in Eqs. (35)- 
(37) does not satisfy Eq. (38), except for special a(x) such as 
a(x) = 1. Indeed, one can show that the coarse-grid operator so 
produced is equal to the right-hand side of Eq. (38) plus some 
additional terms which are due to aliasing effects. This is a 
simple, but lengthy calculation. Here one should be sure to use 
Eqs. (26) and (27) as well as orthogonality relations such as 

N ^ 1 e 2iri jp/N _ 
j=0 

One can achieve a better approximation to the property of 
Eq. (38), i.e., the aliasing terms are far fewer, by a rather 
simple modification of the pseudospectral method. This technique 
is known as the two-thirds rule. It consists of discarding the 
upper third of the frequency spectrum. On a grid with N 
collocation points Eq. (34) is replaced with 


N p = 0 ( mod N ) 
0 otherwise 


(39) 
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ip 

0 


|p| < N/3 

N/3 < I p| < N/2, 


( 40 ) 


PP 


and the interpolation operators become 


R 




l 

N , 
q=- F +1 


e 2iriq(2j - £)/N 


(41) 


and 


P 



q = 




2ir iq ( j - 2Jl)/N 

c • 


(42) 


The price of this modification is that one-third of the collo- 
cation points are wasted. Thus the two-thirds rule version of 
SMG would have to produce a substantial improvement in the 
convergence rate in order to compensate for its reduced 
accuracy. Although no such examples have yet emerged from our 
numerical experiments, the two-thirds rule option may eventually 
prove to be of some use . 
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Chebyshev Series 

The cosine transform matrix C and the 
differentiation matrix D as given in Eqs. (1.45) and 
be left unchanged from Cl]. The analogs to Eqs. (26)- 


“j* " 


Nc, 




[2 

q=0 


1 

c 

q 


cos (ilia, cos (1*2) 


and 


where 


where 


4* " { K 1 cos (t 3 ) cos (^) 


^ Nc^ q=0 q 




c = 

q 


2 

1 


q = 0 or N 
1 < q < N - 1 


c = 

q 


2 

1 


q = 0 or N/2 
1 < q < N/2 - 1 


and 


Rj * ' it. + 


(Q-J_oo Q-4 j.Oo ) i 


j* N S j+2fc 

C £ 


Chebyshev 
(1.47) will 
(30) are 

(43) 

r (44) 

(45) 

(46) 

(47) 

(48) 
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^ + ~ r = 0 (mod N) 

1 . 1 r . 1 , 1. irr>i . ,irr» .nr. .. 

^ + -j cos[(j + jj) -j— Jsin(^-)csc(jj-) otherwise 


( 49 ) 


( N/4 r HO (mod N) 

i- - ^ cos (-j— )+ ^ cos ( (■j + ^)sin(l|)csc(^) otherwise 

(50) 


Equations (43) and (44) represent the "obvious" restriction 
and interpolation operators. Both mey be implemented efficiently 
by Fast Cosine Transforms. Unlike the Fourier series case, 
however, R and P are not adjoint (even aside from a constant 
multiple) unless the boundary conditions happen to be homogeneous 
Dirichlet. A common choice in finite difference multigrid 
algorithms and a natural one in finite element cases is to force 
R to be adjoint to P. Our own computational experience with 
Chebyshev SMG leads us to endorse this strategy. For residual 
transfers, then, we recommend 

R. = N f 2 - 1 coa(i^3)cos(Ii3) , (51) 
Nc. q=0 q N N 


which reduces to 


R 


j* 


Nc j 


(°2j-* + °2-i + 2.) * 


2j+* 


(52) 
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However, the construction of the filtered-coef ficient version of 
the coarse-grid operator via Eqs. (35)-(37) should still be based 
on the restriction formula in Eq. (43). 

IV. AN IMPROVED PRECONDITIONING FOR DIRICHLET PROBLEMS 
Consider the self-adjoint elliptic equation 


d 

9x 


[ a( x,y) 


3u i 

3x J 


+ |y 


3u -I 
3y J 


f 


(53) 


on (-1,1) x (-1,1) with Dirichlet boundary conditions. The 
appropriate pseudospectral approximation employs Chebyshev poly- 
nomials. The collocation points (x.,y ) satisfy 

3 


x . = cos(tt j/N) 

j, £ = 1, 2,. . ., N— 1 . (54) 

y = cos (ir£/N) . 

Let = (N-l)^ denote the total number of degrees of freedom. 

The pseudospectral approximation leads to a discrete set of 
equations like Eq. (6). A detailed description of the matrix L 
representing the Chebyshev discretization of Eq. (53) is given in 
Cl]. 

It is apparent from Eq. (18) that the convergence rate of 
Richardson's iteration on a single grid is governed by the ratio 
of the largest-to-smallest eigenvalues of L. This ratio will be 
referred to below as the single-grid condition number. The 
multigrid condition number, on the other hand, is the ratio of 
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the largest eigenvalue to the smallest high-frequency 

eigenvalue. It controls the smoothing rate (see Eq. (20)). The 

estimates given in [1] for these eigenvalues are X = OfN 1 *), 

max 

X .. = 0(N 2 ), and X . = tt 2 /2 . The implication is that 

mid min 

effective preconditioning is essential for multigrid as well as 
for single-grid iterative schemes. 

Preconditioned Richardson iteration can be expressed as 

v v + u)H~^ (F - Lv) (54) 

where H is a preconditioning matrix . An obvious choice for H is 
a finite difference approximation H FD to the differential 
operator in Eq. (52). In more than one dimension, these finite 
difference approximations are themselves costly to invert. An 
attractive alternative is to use instead an approximate LU- 
decomposition of H FD , i.e., H is taken as the product of a lower- 
triangular matrix 3? and an upper-triangular matrix . In one 
such type of preconditioning (originally proposed by Buleev [6] 
and Oliphant [7] for finite difference discretizations of Eq. 
(53)), denoted by H LU , 3 is identical to the lower-triangular 
portion of H FD and fy/ is chosen so that the two super diagonals of 
LU agree with those of H FD . In [1], a similar decomposition, 
denoted by H RS , was proposed in which the diagonal elements of 3 
were altered from those of H FD to ensure that the row sums of H RS 
and H fd are identical . Incomplete LU-decompositions have been 
used by Wesseling and Sonneveld [8] for multigrid solutions of 
finite difference discretizations. 

Both types of preconditioning can be computed by a simple 

Let (x.,y ) be an interior point of the grid. 

J -v 


recursion. 



Suppose that the finite difference approximation to Eq. (53) at 
this point is given by 


b j,A U j,A-l + d j,£ U j-l,£ + S j,£ U j,£ + f j,£ U j+l,fc + b j,£ U j,£+ 1 f j£ * 

(56) 


The lower-triangular matrix 3? has the non-zero elements 

b. g , d. and e. and the upper-triangular matrix Ql has unit 

diagonal plus the non-zero elements f . and h . where 

3 » £ 3 » £ 


b j,£ b j,£ 


d . = d . „ 

3 / & 3 , £ 


e . „ = 5. . - b . „h. „ . - d. „ f . , „ - ot (b . . f . . , + d . .h . , „ ) 

3 / £ 3»£ 3»£ 3 , £-1 3»£ 3“1»«- 3 * £ 3»£“1 3*£ 3“l/£ 


^j,£ , £^ S j , £ 


(57) 


b j,£ b j / £^ S j # £■ 


The H lu result uses a = 0 and H RS uses a = 1. Straightforward 
modifications are made near the boundaries. 

The eigenvalues of the iteration matrices H“^L corresponding 
to these three types of preconditioning have been computed 
numerically by the QR algorithm [9]. The extreme ones are given 
in Table I. In all cases, the region between 1 and 2.4 is fairly 
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uniformly populated with eigenvalues. In the H LU version there 
are a few (roughly 15%) eigenvalues between X m ^ n and 1; likewise, 
about 20% of the eigenvalues of the H Rg preconditioning fall 
between 2.4 and X . A few of the smaller eigenvalues have 
small imaginary parts. The remaining eigenvalues are real. In 
order to assess the effectiveness of these preconditionings in 
multigrid calculations, one also needs to know the smallest high- 
frequency eigenvalue. The numerical results indicate that this 
is 1.22 for H fd and H LU and 1.45 for H Rg , essentially independent 
of N. The relevant condition numbers are given in Table II. 
Both H l1 j and H Rg require only 0 operations to invert. Thus, 
we reach the striking conclusion that although H Rg is more 
effective for single-grid iterations, H LU is noticeably superior 
in multigrid applications. 


TABLE I 


N 


Extreme Eigenvalues for Preconditioned Chebyshev 
Operator in Two-Dimensions 






X . 
mm 

^max 

^min 

*max 

*min 

*max 

4 

1.000 

1.757 

0.929 

1.717 

1.037 

1.781 

8 

1.000 

2.131 

0.582 

2.273 

1.061 

2.877 

16 

1.000 

2.305 

0.224 

2.603 

1.043 

4.241 

24 

1.000 

2.361 

0.111 

2.737 

1.031 

5.379 
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TABLE II 



Condition 

Number for 

Preconditioned Chebyshev 



Operator in 

Two-Dimensions 



Single- 

-grid 

Multigrid 


N 


H is L 

H 'iu L 

h rs l 

4 

1.85 

1.72 



8 

3.91 

2.71 

1.79 

2.07 

16 

11.62 

4.07 

2.12 

2.92 

24 

24.66 

5.22 

2.26 

3.79 


Beyond N = 24 computations of the complete eigenvalue 

spectra are impractical since the full two-dimensional matrix 

then takes over a million words of storage. The multigrid 

condition numbers and smoothing rates given in Table III are 

based on iterative calculations of the extreme eigenvalues 

of Hr?lL for N = 32 and N = 64 and the empirical formulae, 

LiU 


X = 1.381 N 1 / 8 
max 

X . = 28.37 N -7 / 4 


(58) 


for N > 64. 

The more important of these is the former and it is accurate 
to better than 1% for N > 16. These results suggest that 
0 JlrvfO operations are required for convergence of the SMG 
method based on the H LU preconditioning. This is only slightly 
worse than the best possible result of 0 . The 1- 

parameter smoothing rates are based on a stationary Richardson 
iteration, whereas the 3-parameter smoothing rates are based on 
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nonstationary Richardson iteration employing 3 district 
parameters. Comparing the present smoothing rates with those 
given in [1] for the H RS preconditioning it is apparent that the 
use of the H LU preconditioning reduces the number of SMG 
iterations for N > 32 by at least a factor of 2. 

TABLE III 


N 

Multigrid 

1-Parameter 

3-Parameters 


Condition No. 

Smoothing Rate 

Smoothing Rate 

8 

1.863 

0.301 

0.194 

16 

2.134 

0.362 

0.236 

32 

2.324 

0.398 

0.262 

64 

2.525 

0.433 

0.287 

128 

2.763 

0.469 

0.313 


V. COMPUTATIONAL EXPERIENCE 


The 

global character 

of pseudospectral approximations 

sharply 

distinguishes them 

from local approximations such as 

finite 

difference ones. 

Admittedly, this 

character makes 


spectral methods more complicated to implement, but it also is 
responsible for their superior approximation properties. One 
should expect that somewhat different considerations are 
important in SMG codes than in finite difference ones. Here we 
report on the performance we have obtained with several variants 
of SMG on two-dimensional problems and offer some general advice 
on their use . 

The operation count for a single relaxation sweep is funda- 
mentally different for SMG — OitJV rather than 0£/fO . Thus, 
the standard multigrid accounting device of assessing the cost of 
a coarse-grid relaxation as one-fourth of the work on a grid with 
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half the mesh size is inappropriate. We prefer to make our 
comparisons in terms of actual machine time rather than Brandt's 
work units [10]. This choice has the virtue of including all 
auxiliary effort such as the various interpolations, but it also 
has the disadvantage of depending on the quality of the 
programmer and the computer. 

The specific measure to be used is the equivalent smoothing 
rate, denoted by and defined as follows. In some preliminary 
calculations, the average time t q required for a single fine-grid 
relaxation is determined. For an actual multigrid calculation 
let r*L and be the residuals after the first and last fine-grid 
relaxations, respectively, and let x be the total CPU time. Then 

1 

f- - 1 

o . (59) 

In all the runs reported here, the finest level K = 5. The 
four types of schedules that were examined are described in Table 
IV. In schedules A and B, the problem was first solved on level 
2; then that solution was interpolated to level 3 as the initial 
guess for a multigrid iteration involving levels 2 and 3; then 
the converged level 3 solution was interpolated to level 4 as its 
initial guess, etc. The other two schedules simply began on 
level 5. Most schedules were run in the so-called accommodative 
mode, i.e., the anticipated smoothing rates (e.g., Table III) 
were used in a dynamic determination of when to shift between 
levels. Schedule D used the simple fixed schedule of performing 
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just one sweep through the parameter sequence on a given level 
before interpolating to another one. All runs employed the cor- 
rection scheme CIO] and used random numbers for the initial 
guess. The difference between schedules A and B lies in the 
right-hand side used to define the lower-level problems. In the 
unfiltered version the pointwise values of f(x) were used, but in 
schedule B, the lower-level right-hand sides were obtained by 
applying the appropriate restriction operator to the finest-level 
right-hand side . 


Schedule 


A 

B 

C 

D 


TABLE IV 


Description of Multigrid Schedules 

First Level Control Mode Lower-Level 

Problem 


2 

2 

5 

5 


accommodative 

accommodative 

accommodative 

fixed 


unfiltered 

filtered 

N/A 

N/A 


Periodic Problems 

The test problems have the form of Eq. (53) with the coef- 
ficients 


a(x,y) = b(x,y) = 1 + Ee cos ( B ( X+ V > ■> (60) 


and the exact solution 


u(x,y) = sin(cnr cos x + u/4) sin(air cos y + tt/ 4) - u Q0 , (61) 
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where 


u 00 - \ J o <«> 


(62) 


guarantees that the mean value of u(x,y) vanishes. The source 
term f(x,y) is adjusted accordingly. The parameters of several 
test problems are listed in Table V. The last column relates the 
relaxation parameter actually employed to the optimal parameter 
for the e = 0 problem as given in Eq. (1.37). 


TABLE V 

Parameters of the Periodic Test Problems 


Problem # 

e 

a 

3 

tli/ 

1 

0.00 

1 

1 

1.00 

2 

0.10 

1 

1 

0.75 

3 

0.20 

2 

2 

0.50 


The influence of the coarse-grid operator is indicated in 
Table VI. The filtered coarse-grid operator is defined by Eqs . 
(35) - (37). The unfiltered one replaces Eq. (37) with the 

pointwise values of a(xj). Schedule C was used for all runs. 
The filtered operator presents essentially no improvement. 
Although we find this result puzzling in light of the corre- 
sponding results for Dirichlet problem, we did not pursue it 
further because there are few applications for purely periodic 
boundary conditions . 
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TABLE VI 


Influence of Coarse-Grid Operator on Fourier SMG 


Coarse-Grid 

Operator 


Problem 1 


Equivalent Smoothing Rate 
Problem 2 


Problem 3 


unfiltered 


0.76 


0.78 


filtered 


0.76 


0.78 


0.81 

0.81 


The dependence upon scheduling is given in Table VII. These 
runs used nonstationary Richardson iteration with three distinct 
relaxation parameters as described in Cl]* The filtered coarse- 
grid operators were also employed. The most striking result is 
the distinct superiority of schedule B on problem 1. The expla- 
nation for this behaviour lies in the very special relationship 
that exists for the constant-coefficient problem between the 
interpolation operators R k and P^ and the operators L^: The 
eigenfunctions of are a subset (in fact precisely the low- 
frequency subset) of the eigenfunctions of L^. The prolongation 
operator P^ leaves the eigenfunctions of unchanged. Thus, 
this interpolation introduces no spurious high-frequency com- 
ponents. A similar relationship holds for the restriction 
operator. 
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TABLE VII 


Influence of Scheduling on Fourier SMG 


Schedule 

A 

B 

C 

D 


Equivalent Smoothing Rate 
Problem 1 Problem 2 Problem 3 


0.56 

0.50 

0.61 

0.70 


0.70 

0.70 

0.70 

0.76 


0.82 

0.82 

0.80 

0.83 


In order to get the full benefit of this property, however, 
the lower-level problems (used for obtaining initial guesses on 
the higher-level problems) must have alias-free right-hand 
sides. Consider how the simple model problem described by Eq. 
(1) behaves in transform space for an elementary multigrid scheme 
which uses only levels 2 and 3. The (transformed) level 3 
equations are 

p 2 u p = fp p = -4, -3,.., 3. (63) 


Suppose that a fully-converged level 2 solution is used for the 
initial guess on level 3. If the level 2 problem is defined by 
the restricted values of f(Xj) (as in schedule B) , then its 
equations are 



p = -2, -1,0,1. (64) 


The resulting interpolated initial guess for the level 3 problem 

A A 

will have u p for p = -2, -1,0,1 precisely correct and u p for p = 
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-4, -3, 2,3 will be 0. Thus, the only errors in the level 3 
solution will be in the high-frequency modes and there will be no 
need to make any coarse-grid corrections, i.e., no true multi- 
gridding will occur. On the other hand, suppose that the level 2 
problem is defined by the pointwise values of f(Xj) (as in 
schedule A). Then the level 2 equations will be 



p = -2, -1,0,1 , (65) 


where 



( 66 ) 


The last term in each of these equations is, of course, the alias 
of its preceding term. When this converged level 2 solution is 
interpolated to level 3 for use as the fine-level's initial 
guess, there will be errors in the low-frequency modes. Hence 
coarse-grid corrections will have to be made. Consequently, 
schedule A consumes more computer time than schedule B. 

The superiority of schedule B does not extend to non-trivial 
cases, as represented here by problems 2 and 3. Since the 
eigenfunctions of the discrete operators are no longer simple 
trigonometric functions, they are not preserved by the interpo- 
lation operators. 
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Schedule D is clearly marked as inferior to schedule C. 
Although schedule C operates in the accommodative mode, in nearly 
every case there are 6 relaxations on a level before restriction 
to a coarser level occurs. Schedule D uses only half as many 
relaxations before restriction occurs. (Recall that 

nonstationary Richardson iteration with 3 parameters is employed, 
so that the number of relaxation is necessarily a multiple of 
3.) Part of the increased efficiency of schedule C arises from 
the less frequent use of interpolations. 

A Fourier SMG program has some subtleties that deserve 
mention. They are connected with the zero eigenvalue of the 
discrete operator that arise from the p = 0 and p = - ^ diagonal 
entries of the Fourier differentiation operator (see Eq. (34)t. 
The associated mean-value and highest-frequency eigenfunctions 
are undamped by the iterative scheme. For the constant-coef- 
ficient case, a sufficient precaution is to filter these com- 
ponents out of the right-hand side and the initial guess. In 
variable-coefficient problems one must ensure that none of the 
highest-frequency component enters the solution during a 
relaxation. 

Dirichlet Problems 

The test problems are specified by 

a(x,y) = b(x,y) = 1 + e e cos ( ( x+y ) ) < 67 ) 

u(x,y) = sin(cnrx + tt/ 4) sin(airy + tr/4). (68) 
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The parameters of the test problems are given in Table VIII. 
Problem 1 not only has constant coefficients but it is also so 
well-resolved by the Chebyshev pseudospectral method that its 
discretization error on level 5 is well below the round-off error 
of the CDC Cyber-175 (14 digits). Problem 3 is at the other 
extreme. The coefficients of the equation oscillate so rapidly 
that the level 5 grid cannot resolve them. Instead the converged 
solution of the level 5 collocation equations has an error of 
order 1. This case is included as a test of whether the 

Chebyshev SMG method is robust enough to converge on such a 
problem. 


TABLE VIII 


Parameters of the 


Problem # e 

1 0.00 

2 0.20 

3 1.00 


Dirichlet Test Problems 

a 3 oj/oj* 

1 1 1.00 

2 2 1.00 

5 10 0.90 


The difference between the two choices of the coarse-grid 
operators is shown in Table IX. Both versions are identical on 
problem 1. The unfiltered coarse-grid operators produce a slow 
method on problem 2. The extra work occurs on the coarser levels 
where the smoothing is less effective. On problem 3 the 
unfiltered coarse-grid operators lead to a divergent method. We 
find it curious that the filtered coarse-grid operator failed to 
produce a similar improvement in Fourier SMG. Perhaps the 
difference lies in the use of preconditioning for Chebyshev 
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SMG. (The filtered coefficients are, of course, used in the 
finite difference preconditioning as well.) Note that the 
variable coefficients of problem 3 are extremely oscillatory. 
Our unpreconditioned Fourier SMG method cannot handle anything 
remotely as difficult. 


TABLE IX 


Influence of Coarse-Grid Operator on Chebyshev SMG 


Coarse-Grid 

Operator 


Problem 


Equivalent Smoothing Rate 
1 Problem 2 Problem 3 


unfiltered 


0.62 


0.66 


divergent 


filtered 


0.62 


0.62 


0.76 


The scheduling dependence is given in Table X. The filtered 
coarse-grid operator was used along with nonstationary Richardson 
iteration. The same trends are apparent here as for the periodic 
case, except that schedule B is now only slightly better than 
schedule A on problem 1. In the Chebyshev method the 
eigenfunctions of the discrete, constant-coefficient operator are 
not preserved by the interpolation procedures. As before 
schedule D nearly always produces 6 relaxations prior to 
restriction. We again see that it is not a good strategy to 
relax the minimum number of times before restricting. 
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TABLE X 


Influence of Scheduling on Chebyshev SMG 
Equivalent Smoothing Rate 

Schedule Problem 1 Problem 2 Problem 3 


A 0.26 

B 0.25 

C 0.51 

D 0.60 


0.58 

0.77 

0.57 

0.75 

0.59 

0.70 

0.67 

0.74 


At the end of the third section, we recommended the 
restriction operator defined by Eq. (51) for the residual 
transfer. Two alternatives have been tested. The "obvious" 
restriction operator of Eq. (43) fails miserably. It doesn't 
even work for the constant-coefficient problem. In the 
accommodative mode the algorithm rapidly settles into a "limit 
cycle" involving levels 2 and 3: it alternates between these two 
levels, always arriving on either level with the solution in the 
same state as at the start of its last visit. The second 
alternative is to "homogenize" the restrictions and prolongations 
by forcing the boundary values of the corrections to be zero both 
before and after the interpolation. Although this bizarre choice 
was made by accident, it actually works. However, since it has 
uniformly been slightly less effective than the adjoint choice, 
there is no good reason to resort to it. 


VI . PERSPECTIVE 

Together with [1] this paper represents a comprehensive 
description of the fundamentals of the spectral multigrid method 
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for solving pseudospectral discretizations of self-adjoint 
elliptic equations. The key elements of the coarse-grid 
operators and the interpolation operators have been described in 
detail and their efficacy demonstrated in the numerical 
examples. No doubt the spectral multigrid method will prove as 
amenable to improvements in the relaxation scheme and precon- 
ditioning aspects as have finite difference and finite element 
multigrid methods. The next big step is the demonstration of the 
power of this method on a difficult, nonlinear problem of 
engineering interest. This has already been achieved and is 
reported in [11]. 
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